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On ■ We study nonequilibrium properties of a one-dimensional lattice 

>~>' Hamiltonian with quartic interactions in strong thermal gradients. 

Q ■ Nonequilibrium temperature profiles, T{x), are found to develop signif- 

c^ : 

i-C ■ icant curvature and boundary jumps. From the determination of the 



bulk thermal conductivity, we develop a quantitative description of T(x) 
including the jumps. 
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The description of transport in physical systems is usually relegated to the near equilib- 
rium regime, where Green-Kubo theory can be used. When one strays from this into systems 
far from equilibrium, far less is known about the statistical mechanics, or the behavior of 
transport coefficients |I|]. While the physics of non-equilibrium systems is certainly of broad 
interest, many basic problems have yet to be fully understood. Of particular interest are the 
non-equilibrium stationary states and the nature of the statistical mechanics which char- 
acterize it. One approach to understanding this problems is to place systems in thermal 
gradients and examine the long time behavior of observables. In one dimension, this has 
been done in a variety of problems where both finite and divergent transport coefficients were 
measured. Typically the transport coefficients are found to diverge when the Hamiltonian 
conserves total momentum. This is typical of systems where the interactions depend only 
upon differences Xi — Xj, such as the FPU and Toda chains [0-0. When an 'on site' poten- 
tial, V{xi), is also present, the coherent propagation of long wavelength modes is suppressed 
resulting in finite conductivity Q. This is the case in the Frenkel-Kontorova, ding-a-ling, 
Lorenz and other models H-B[. Bulk behavior is also known in higher dimensions as well 



In this letter, we investigate the non-equilibrium steady state properties and thermal 
transport of lattice Hamiltonians with quartic interactions in 1 spatial dimension. Our 
system is the discrete version of 0^ scalar field theory, a proto-typical field theory that has 
broad application. In contrast to systems such as the FPU jS model P,3JT0(], which have 
divergent thermal conductivities, we find a well defined bulk limit for our non-equilibrium 
results. This we attribute to the on-site nature of the 0^ interaction. Certain properties 
of this theory have been studied in the past, which include the Hamiltonian dynamics and 



phase transitions, as well as the ergodic properties |Tl|. However, the thermal conductivity 
has not been determined, and further, there is yet no understanding of the role of boundary 
jumps and its inter-relation with the shape of the non-equilibrium thermal profile. We 
present a quantitative analysis of this effect which describes the behavior near and far from 
thermal equilibrium. A full description of the temperature profiles far from equilibrium is 



shown to require a description of the boundary jumps. 
Our model system is the 1-dimensional Hamihonian 
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where /?^, ?7i and x ^i-re parameters. This model becomes identical to the FPU /? model if 
we substitute [qi+i — gj)^ for g^^ in the interaction. It is convenient to perform the rescalings 
Qi = Qi{x/P)^ Pi = PiiV^X^/P) and t = t^y/m/x)- In the dimensionless variables qi,Pi,t, 
the Hamiltonian is 



^ i=\ 



Pi + (?*+i - li) + 2^i 



(2) 



where B. = H'f3'^/x^- To study the statistical mechanics of H in strong thermal gradients, 
we evolve Hamilton's equations of motion together with thermal boundary conditions on 
the two ends. Observables are allowed to converge to steady-state values, which are then 
analyzed. The dynamics is solved on a spatial grid using either fifth and sixth order Runge- 



Kutta or leap-frog algorithms |]T2|. We impose time-reversal invariant, thermal boundary 
conditions on the equations of motion at sites i = 1 and i = L. Specifically, the endpoint 
equations are augmented to include dynamical thermal constraints using the robust methods 
of [|I^,|I^, and become: 



(3) 



qi = Pi, qi = PL, 

P^ = -If - t^?Pi - t^2P?, PL = -If - ^wIpl - ^fw^pl, 

w, = ai(p2/r0 - 1), ws = as{pl/T^ - 1), 

W2 = a2ipt/T^ - M) w, = aM/n - ^pD- 

The auxiliary variables Wk dynamically implement the statistical boundary conditions con- 
sistent with the boundary temperatures T° and T° [|13|. We will use the optimal choices 
for the couplings a,, which in this case are Oj = 1 [|14|. ( We have checked that the results 
presented here do not depend on this particular choice of the thermostats Wk, the values 
of Oj, nor on the number of sites at each end that are thermostatted. We have also exam- 
ined thermostatting from 1 to 3 sites on each end; we return to this point when we discuss 



boundary jumps.) With this control of the endpoint temperatures, we are able to simulate 
the dynamics of H with the boundary thermostat temperatures (T°, T°) fixed. Apart from 
the endpoints, the system evolves according to the dynamics dictated by the Hamiltonian 
(0). In the simulations, we use between 10^ and 10^ with time steps of dt from 0.1 to 0.001. 
The lattice size was varied from L = 10 to 8000. 

In order to understand the transport properties, it is useful to start with the stress-energy 
tensor T^'^ . From the continuity equation, 

.T^-^= E + VJ = 0, (4) 



dxf^ dt 

where E = T^^ is the energy density. For thermal gradients in the x— direction, we can 
identify J = T^^{xi,t) = —pi{qi+i — Qi) as the heat flux. 

T° = T^ : By setting T° = Tj^ = T, we verify that the canonical ensemble is realized at all 
sites. This is done in several ways. By following the trajectories Pi{t) and histogramming 
the values at each time step, one reconstructs the entire momentum distribution at any site, 
which are found to converge to f{pi) oc exp[— pf/2T]. We also verify the virial theorem for 
the kinetic and potential energies. 

In thermal equilibrium, we can use the Green-Kubo formula to compute the thermal 
conductivity k{T) in the linear-response regime: 

1 r°^ 
k{T) = —- dtJ2{J{x,,t)J{x,,,0))^^. (5) 

^^-^ •'^ k,k' 
The autocorrelation function is evaluated in the canonical ensemble, T° = T2, and the 
number of points used in the sum is A^ < L since we omit points near the boundaries. The 
formula is expected to hold within the linear regime, at least in higher dimensions. In 1-d, 
the integrand in (^ has been argued to develop a long time tail ~ t^^/^ on time scales on 
the order of a few ten times the mean free time, leading to the divergence of (^) [|T^ . In our 
case, the presence of the on site interactions serve to destroy this long-time tail at sufficiently 
large times. We plot the absolute value of the integrand of (^) for a temperature T = 1/10 
in Fig. 1 (top). At this temperature the mean free time roughly 50, so that long-time tails 



should be present on times ~ 10^. The t~^/^ behavior is indicated by the dashed hne. The 
time integral is given in 1 (bottom), showing that the integral (|^) is finite. These linear 
response predictions will be seen to agree with the direct measurements we discuss below. 
We see that on the time scales of t ~ 50, there is some agreement with the long-time tail 
predictions. However, on longer times the divergent behavior is not present, and the linear 
response results converge to a well defined limit. In Fig. 2, we compile the Green-Kubo 
measurements of k, plotted as a function of T. In the context of the FPU model, it has 
recently been shown that the long-time tails of the autocorrelation function behave as t~^/^ 
leading to a t^^^ divergence, based on mode-coupling theory [|TB|. While we do not have a 
divergence in our integrals, the behavior of the integrated autocorrelation function in the 
transient regime is much closer to t^/^ at all measured temperatures. 

T^ > T° : We can verify the linear response calculations through a direct measurement of 
the heat fiow. The thermal conductivity n is defined using Fourier's law {J)ne = —i^dxT{x), 
where {J)ne is the averaged heat fiux. The gradient is evaluated by taking |T° — T°| suf- 
ficiently small so that the temperature profile is linear. We then vary the temperature 
difference |T^ — T°| around the same average temperature to verify that J is proportional to 
VT, and extract n{T). T{x) is the local temperature defined through an ideal gas thermome- 
ter, by T(x) = (p^(x))jvB, where p{x) is the momentum at site x. Here (■ ■ ■)^rB indicates the 
ensemble average over the non-equilibrium steady state. To obtain the transport properties, 
each simulation is run long enough for observables such as J, the energy density as well as 
distribution functions to converge. In Fig. 2 we show the results of both direct measure- 
ments and Green-Kubo predictions for n. We find that both agree over several orders of 
magnitude, and that k{T) has the power law temperature dependence 

«:(T) = A, 7 = 1.35(2), A = 2.83(4), (6) 



reminiscent of the behavior of lattice phonons at high temperature IT^. We have also 



verified that a sensible bulk behavior exists, as shown in Fig. 3; the thermal conductivity 
is independent of L when it is larger than the mean free path, which, on the lattice, is of 



order of the conductivity (see below). The dashed hnes in Fig. 3 correspond to the values 
of Eq. d^) at that temperature. 

T^ ^ T° : By controlling T^ and T° we can begin to explore the non-equilibrium steady 
state. For T° ~ T^, we expect to be in the linear regime, with a linear temperature profile, 
and transport given by Fourier's law. This has been readily verified. As the difference be- 
tween endpoint temperatures increases, two characteristics emerge: the temperature profile 
develops curvature (Fig. 4(a) where T^/T^ = 0.05), and there are substantial temperature 
jumps near the boundaries (Fig. 4(c)). 

When the temperature varies substantially in the system, one cause for the non-linearity 
is the temperature dependence of the thermal conductivity. If we assume that the tempera- 
ture dependence of the thermal conductivity is the dominant source of the non-linearity of 
the thermal profile, T(x) can be obtained by integrating Fourier's law as 



T(x) = Tr_ 



-'l)'1f 
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, 7 7^ 1. (7) 



The integration is simple since the energy flow {J)ne is independent of x due to current 
conservation. This result is distinct from temperature profiles discussed previously as in 
Ref. [0. In this equation, the temperatures Tc^h denote the temperatures at the boundaries 
obtained by extrapolating the temperature profile. These are in general different from the 
endpoint temperatures T\. (In Fig. 4(c), T° = 0.05 and Tc = 0.096.) The temperature 
profile is a function only of x/L so that it has a smooth continuum limit. The ratio of Eq. 
(^) to the measured local temperature is shown in Fig. 4 (b). Aside from the endpoints, 
one can see that the agreement is quite good. We have found that this formula works well 
in higher dimensions as well (rf = 2, 3) [p^ . 

It is clear that while the formula for T{x) describes the shape of our observed non- 
equilibrium profiles, it depends on quantities which are not determined from our input 
parameters, namely the extrapolated temperatures. Hence a complete description of T{x) 
will require the understanding of how the boundary jump Tc — T° (and similarly for T^) 
depends on the parameters of the system. We will do this below. 



From Eq. (|^), we find the heat flux in the non-equilibrium steady state to be 
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We can expand this around the average temperature Tav = (Th + Tc)/2, to find: 
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Here {J)%e i^ identified as the constant (temperature) gradient limit of the heat flux, which 
actually provides a good approximation to the heat flux even when there is curvature in the 
temperature profile, providing AT /Tav ^ 1- For AT /Tav ~ 1, the notion of local equilibrium 
becomes questionable, and this description of the heat flux starts to break down [Q. 

As we mentioned previously, a common feature to non-equilibrium steady states is the 



appearance of boundary jumps or slips &|2,19]. These are not artifacts of the simulation, but 



are true physical effects [19,20|. For systems in thermal gradients, the boundary temperature 
can be different from the temperature of the system near the edges. Similar behavior can 
also appear in fluids which are sheared [|2T]]. In this case there is slippage between the system 
and the wall. Such effects are well known sources of error in the experimental measurement 
of the transport coefficients in fluids [^. In our lattice model, elementary kinetic theory 
predicts a thermal conductivity of k = CyCgl for the system, where Cy is the heat capacity 
per volume, Cg the sound velocity and I the mean free path. The formula applies to the 
basic excitations, the "phonons", of the system. In our model, Cv,Cs ~ 1. Therefore, the 
mean free path is expected to be of order of the conductivity itself. This allows us to clarify 
the nature of the temperature jumps at the boundaries. The jump, arising due to the finite 
mean free path, satisfies a relation 



T.-T.° = ,^ 



dT 

dx 



(10) 



boundary 



boundary 

when I <^ L |]19|. Here dT /dn is the normal derivative at the surface, so that the + (— ) sign 
corresponds to the lower (upper) edge of the chain. We will consider the lower temperature 
end (the + sign), but identical results hold for the high temperature side if one changes the 



sign accordingly, rj is expected to behave like the mean free path and hence the conductivity, 
roughly speaking. If we let rj = aK for some constant a, we can use Eq. (|^) to find the 
approximate behavior of the boundary jumps: 

T,-T!^c:^a{J)NE (11) 

rrpO rpQ\CX,K[l av) 
~ {J-h - J-c) ^ \ • 

We have determined tj by fitting Eq. (^ to profiles T[x) obtained from various non- 
equilibrium steady states, using the power 7 obtained from Eq. (^. The behavior of rj 
as a function of T from the non-linear thermal profiles is plotted in Fig. 5 (top), which 
can be fitted to 77 = (6.1 ± 0.5) x 7"-i-5±o-i^ consistent with the above argument. We can 
also use (|lT]) to study the boundary jumps directly. In Fig. 5 (bottom), we see that Eq. 



(pAl) is readily verified in the data, which includes thermal profiles both near and far from 
equilibrium. A simple fit (dashes) provides the coefficient a = 2.6(1), fully consistent with 
r]{T) = aK,{T). We have analyzed these effects in higher dimensions as well and find similar 
results ig. 

We now have the following behavior: Thermal boundary conditions will result in a 
given amount of heat fiux, which in turn determines the magnitude of the boundary jumps 
through Eq. (11). Once the boundary jumps are defined, the thermal profile, Eq. (7), is 
determined from the endpoint temperatures T°. We have examined the dependence of the 
results we present here on the types of thermal boundary conditions we have used. This 
includes variations not only of the strengths of the interactions, Oj, but also of the number 
of thermostatted sites on each end. We have observed that the coefficient a (and hence r]) 
is not a physical property of the Hamiltonian, but can depend on the thermostats. Different 
thermostats can produce different boundary jumps for the same boundary temperatures 
T\. Although the new jumps provide new extrapolated temperatures Tc^h, the temperature 
profile is still predicted by Eq. (7). Hence the formulas we have presented are still valid 
in spite of the (slight) differences in boundary jumps. These differences do not lead to any 
modification of the transport coefficients, and are only differences in the shape T{x) due to 
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the variation of the boundary jumps which seem to reflect the efficiency at which heat can 
be transferee! to the system through the boundary in the far from equihbrium regime. 

We have constructed the non-equihbrium steady states of an anharmonic chain with 
quartic interactions, by placing it in strong thermal gradients, and obtained well defined 
bulk transport properties. These results are also valid for classical 0^ lattice field theory. 
We determined the temperature dependence of the thermal conductivity, and then derived 
the non-equilibrium thermal profile T{x), which agrees well with the observed behavior in 
the non-equilibrium steady state. A simple expression for the non-equilibrium heat flux is 
also found to agree well with direct measurements. Using these results, as well as arguments 
from kinetic theory, we are able to quantify the temperature jumps at the edges of our 
system, which are endemic to both non-equilibrium simulations and experiment. We hope 
to develop this overall picture of the non-equilibrium steady state towards understanding 
the dynamics of phase transitions under non-equilibrium conditions, as well as the concepts 
of local equilibrium. 
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FIG. 1. (Top) Absolute value of the autocorrelation function 

C(.t) = J2k,k'{J{xk,t)J{xk',0))EQ/NT'^, for L = 100, up to a time t = 2 x 10''. The long-time 

tail divergence would have the behavior shown by the dashed line. (Bottom) Green-Kubo integral 

up to time t, showing convergence of k. The dashed line is the anticipated behavior if the long-time 

tail divergence were present. One can see that on the time scales t ~ 10^ (a few ten times the 

mean free time), similar behavior can be seen, although it vanishes for larger times. 
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FIG. 2. Thermal conductivity k obtained from direct (□) and Green-Kubo (x) measurements 
for various lattice sizes L, and the power law fit k{T) = 2.83(4)/T"'^''^^'^^^ (dashes). Green-Kubo 
integrals readily converge to values consistent with direct measurements. 
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FIG. 3. Size (L) dependence of the thermal conductivity indicating bulk behavior, for tem- 
peratures (upper to lower) T = 1/10 (+), 1/4 (□), 1/2 (x) and 1 (O). The dashed lines are the 
predictions from Eq. (6) for that temperature. 
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FIG. 4. (a) Non-equilibrium therm al profile T{x) as a function of position (solid), compared 
to theoretical predictions of Eq. (7) (dashes) which lie on the solid curve. The temperature was 
sampled 10® times every At = 0.25. Endpoints were thermalized at {T^,T^) = (0.05, 1). (b) The 
ratio of the predicted profile, Eq. (7), to that measured in (a). Aside from the boundary jumps, 
the agreement is within a few percent, (c) Blowup of the temperature jumps at the low end of (a). 
The simulated temperature T^ can differ significantly from the extrapolated value Tg. 



10^ 



10- 




0.4 




-0.10 



0.00 



-0.05 

<J>]Vjg 

FIG. 5. (Top) Temperature dependence of the jump parameter ij and its approximate power 
law behavior (dashes). (Bottom) Verification of Eq. (11). We plot the measured boundary jumps 
as a function of the non-equilibrium heat flux. The slope provides the coefficient a = 2.6(1) which 
relates ry to the conductivity k via rj = ok. 
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